Probing molecular free energy landscapes by periodic loading 
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Single molecule pulling experiments provide information about interactions in biomolecules that 
cannot be obtained by any other method. However, the reconstruction of the molecule's free energy 
profile from the experimental data is still a challenge, in particular for the unstable barrier regions. 
We propose a new method for obtaining the full profile by introducing a periodic ramp and using 
Jarzynski's identity for obtaining equilibrium quantities from non-equilibrium data. Our simulated 
experiments show that this method delivers significant more accurate data than previous methods, 
under the constraint of equal experimental effort. 
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Introduction - A key feature of biological systems is 
the high degree of self-organization of polymers, pro- 
teins, and other macromolecules, and their interaction 
with smaller components such as energy providers or 
messenger molecules 0. These processes are ultimately 
driven by specific and tunable molecular interactions. 
Their detailed knowledge is thus a prerequisite for the 
understanding how biological systems work on molecu- 
lar and higher levels. Recent developments of highly 
sensitive force probes such as atomic force microscopy 
(AFM) HQ, optical and magnetic tweezers HH@, and 
biomembrane force probes U@] make it possible to probe 
the molecular interactions of individual biomolecules by 
their response to mechanical stress (see jj| for re- 
views). The systems studied by single molecule pulling 
experiments can be divided in two groups: in rupture 
experiments, receptor and ligand molecules are attached 
to a substrate and a transducer, respectively, often via 
chemical linkers. After allowing receptors and ligands to 
bind, the transducer, e.g., an AFM cantilever, is pulled 
away, which causes the receptor-ligand pairs to rupture. 
The maximum force the molecule can withstand has been 
measured in this way for biotin and streptavidin 0, Q , 
and many other receptor-ligand pairs . Secondly, un- 
and refolding experiments probe the elastic properties 
of an individual biomolecule. The molecule is attached 
between a substrate and a transducer, again via chem- 
ical linkers. Force-extension relations are obtained by 
measuring the force as a function of the position of the 
transducer. In this way one may explore regions of the 
free energy landscape of the biomolecule far away from 
thermal folding pathwa ys. Investigated systems include 
DNA [HElEi, RNA [il El, polysaccharides [HI , the 
muscle protein tit in [18| . the membrane protein bacteri- 
orhodopsin (BR) |l9j. and many other proteins |20j . 

Figure ^ shows a typical setup of single molecule 
pulling experiments. The molecule is attached between 
the substrate surface and the cantilever tip. The position 
of the cantilever x(t) is moved according to a prescribed 
experimental protocol. The extension of the molecule is 
described by a suitable reaction coordinate z given by 
the position of the cantilever tip. For fixed extension z 
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Figure 1: Schematic view of the experimental setup and a 
generic free energy potential G(z). The first minimum repre- 
sents the folded state, whereas the second shallow minimum 
represents the unfolded state of the biopolymer. The coordi- 
nate x denotes the position of the cantilever and z the posi- 
tion of the cantilever tip to which one end of the biopolymer 
is attached. 



and time t, the energy of the molecule perturbed by the 
cantilever spring is given by the time-dependent Hamil- 
tonian 



H(z,t) = G(z)+V (z,x(t)) 



G{z) + ^{x{t)-z) 2 



(1) 

where G(z) is the free energy profile of the unperturbed 
molecule. The second term describes the external force 
acting on the molecule in terms of a harmonic potential 
with effective spring constant k. Since the molecule is 
coupled to a heat bath at temperature T the time evolu- 
tion of z is stochastic. 

Traditionally, the cantilever is moved according to a 
linear ramp, 



"(t) = Xq + vt 



(2) 



with offset xq at t — and constant velocity v, and the 
force F(t,v) acting on the cantilever is recorded. The 
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challenge is to recover from these data the unperturbed 
molecule's free energy profile G{z) containing the desired 
information about the molecular interactions. Evans 
and Ritchie first pointed out that the rupture force of 
receptor- ligand pairs depends on the loading rate v |2lj| . 
Thus, by combining the data F(t,v) for a broad spec- 
trum of loading rates v, referred to as dynamic force 
spectroscopy, important features of G{z) can be deter- 
mined such as the distance between the minimum and 
maximum of an energy barrier for rupture |2lj | . Heymann 
and Grubmuller refined this technique and obtained the 
heights and positions of the maxima of a molecular force 
profile d z G(z) with high spatial resolution j^. On the 
other hand, traditional experimental protocols like the 
linear ramp J5J) still entail certain drawbacks. First of 
all, the thermodynamically unstable (concave) barrier re- 
gions of G(z), determined by specific molecular interac- 
tions and therefore of particular interest, are poorly sam- 
pled due to snapping motion and thus hard to determine 

Periodic loading - In an effort to improve the quality of 
data obtained by single molecule pulling experiments, in 
this work we propose a new method for obtaining the full 
free energy profile G{z) by introducing a periodic ramp, 

x(t) = xq + asin(u)t) , (3) 

with given offset Xq, amplitude a, and frequency u>. Fig- 
ure|21shows that periodic loading delivers significant more 
accurate data for the sample free energy profile G(z) than 
the linear ramp (0) , under the constraint of equal experi- 
mental effort. The improvement of the quality of data in 
the important barrier region of G{z) is striking. The bet- 
ter performance of the periodic loading method as com- 
pared to linear loading is mainly due to the fact that pe- 
riodic loading ensures that the barrier region is traversed 
often and from both sides. The quality of our recon- 
struction moreover depends crucially on the fact that we 
sample the barrier region under non-equilibrium condi- 
tions taking advantage of Jarzynski's relation to recover 
the equilibrium profile [23L |24| . Driving the system out 
of equilibrium is important since under quasi-static con- 
ditions an efficient sampling of the barrier region (where 
G{z) 3> ksT) is inhibited by the equilibrium Boltzmann 
factor exp[— G(z)/(fcsT)] <C 1. For periodic loading, 
the freedom to choose the frequency u> large enough al- 
lows one to probe the same region under non-equilibrium 
conditions, thus overriding the exponential punishment 
by the equilibrium Boltzmann factor. The optimal fre- 
quency arises from balancing competing effects as quan- 
tified in a case study below. This frequency should not 
be too large in order to enable the system to follow the 
external drive. Moreover, the Jarzynski procedure con- 
verges the slower the further one is away from equilibrium 
|25| . For too small u>, on the other hand, one does not 
generate enough crossings under the constraint of a finite 
total measuring time. 
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Figure 2: Comparison of reconstructed free energy profiles 
G{z) by using periodic (0) vs linear loading (+), generated 
by JJy. The solid line is the original free energy profile G{z). 
For both methods, 10 trajectories of 75 ms length, a spring 
constant k = 11.6 pN/nm, and a diffusion constant D — 
10 _7 cm 2 /s were used [see Q and Q]. For periodic loading 
we used the optimal frequency uj* — 1.2 x 10 3 1/s (determined 
in Fig.|3J, the preloading offset xo = 6 nm, and the amplitude 
a = 5 nm in For linear loading we used the optimal 

velocity v* = 200 nm/s (determined in Fig.[3J. The inset 
compares the relative deviation of the reconstructed data from 
the original profile G(z) corresponding to the zero line. 

Simulated experiments - We have tested our proposal 
of periodic loading with simulated experiments and com- 
pared it with the traditional method of linear loading. To 
this end we have chosen a generic free energy profile G(z) 
for the unfolding of tertiary structures of biopolymers 
such as the membrane protein BR Our sample free 
energy profile has two separated minima, one of which is 
narrow and deep representing the folded state and one 
of which is shallow representing the unfolded state, see 
Fig-E For BR, a rich structure of unfolding transitions 
under force was found |l9|. Single force peaks in the un- 
folding spectra could be allocated to specific changes in 
molecular configuration. Most of the force peaks scatter 
between 25 up to 100 pN. With a typical length scale of 
several nanometers this yields an energy barrier of about 
20 UbT at room temperature. We focus on one of such 
transitions and choose a barrier of 2 nm length and 23 
ksT height, leading to a typical transition force of about 
50 pN in our simulated pulling experiments. 

For comparing the periodic with the linear ramp, we 
simulated both kinds of protocols using typical parame- 
ters as given in the caption of Fig. [2 We have generated 
an ensemble of trajectories ((t) of the reaction coordinate 
z by discretizing the Langevin equation 

dz/dt = -^ 1 dH(z,t)/dz + V2D^(t) (4) 

with the Hamiltonian H(z,t) from QJ. The Gaussian 
random force £(t) has zero mean and short-ranged tern- 
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poral correlations (£(£)£(£')) = S(t — t'). The diffusion 
constant D is related to the friction coefficient 7 by the 
Einstein relation D = kgT/^. In our simulations, the 
length of a time step is limited by the condition that 
a spatial step should be small compared to the typical 
length scale set by the free energy profile. In addition, 
the recording rate of the data ((t) should be much smaller 
than an inverse time step but large enough to resolve the 
cantilever motion. 

The reconstructed free energy profile for both proto- 
cols is shown in Fig-EI The overall quality of the data 
obtained by the periodic ramp is far better than the lin- 
ear ramp, especially in the barrier region where the data 
obtained by the linear ramp underestimates the barrier 
height by several k B T. In order to ensure an unbiased 
comparison, we have chosen the same total number of 
trajectories and the same total measuring time for both 
methods. In practice, the measuring time of the trajecto- 
ries, once prepared, is not a limiting factor. The periodic 
ramp therefore allows measuring as many transitions as 
necessary to collect the sufficient amount of data. 

Reconstruction using Jarzynski's relation - Pulling pro- 
tocols both in real experiments and in simulations used 
here typically generate non-equilibrium data from which 
one has to recover an equilibrium property like G(z). 
Furthermore, as outlined above, the method proposed 
here purposely takes advantage of the non-equilibrium 
conditions generated by a large enough optimal driving 
frequency uj*, to be determined below. The difficulty 
to recover equilibrium properties from non-equilibrium 



Using this expression, we have generated the data shown 
in Fig.EJ 

Optimization with respect to frequency /velocity - The 
quality of the reconstructed free energy profile depends 
crucially on the driving frequency u> for the periodic ramp 
and the velocity v for the linear ramp, respectively. To 
quantify this observation, we calculate the mean square 
error a 2 = ([G(z) — G(z)] 2 ) Zi where the average is taken 
over discrete values z, in the z-interval under consider- 
ation. For clarity, we denote by G{z) the reconstructed 
free energy profile based on JBJ . For a better reconstruc- 
tion quality, a 2 is smaller. The z-interval is chosen from 
the first to the second minimum. The total measuring 
time and the number of trajectories are the same as be- 
fore. 

Figure shows the mean square error a 2 of the re- 
constructed free energy and some characteristic confi- 
dence intervals (error bars) for both protocols. For peri- 



data may be resolved by using a recent advance in non- 
equilibrium statistical mechanics due to Jarzynski |23| . 
according to which the equilibrium profile G(z) can be in- 
ferred by suitably averaging non-equilibrium trajectories 
((t) of the reaction coordinate z [24| . This method was 
verified by stretching RNA reversibly and irreversibly be- 
tween two conformations indeed 26] . 

Jarzynski's relation, in general, states that the free en- 
ergy difference AG between two equilibrium states can 
be extracted from averaging the work W required to 
drive the system from one state to the other according 
to e _/3AG = (eT m ) with = l/k B T the inverse tem- 
perature [2i|. The generalization from two states to a 
z-resolved free energy profi le G(z) perturbed by a har- 
monic spring (JJJ reads |24j 

e -pmz,t)-ao] = (s [z _ m e -pw(t) \ (5) 

The average (...) is over infinitely many realizations 
< t' < t, of the stochastic trajectory of the biopoly- 
mer's end position, starting in equilibrium at t' = and 
ending at the given position z at if = t, as enforced by 
the delta function. The external work is a functional of 
C(i') and given by W{t) = f* db' d T H(((t'),T)\ T = t , . The 
constant G = -k B T1n(J™ dz' e -^( 2 '<°) /A T ) is the free 
energy in the initial state at t = 0, where the thermal 
wavelength Xt = h/ \Z2irmkBT serves for normalization. 

Summing up the normalized distributions obtained 
from JSJ) at each time slice with the method of weighted 
histograms |2j| yields the reconstruction formula for the 
unperturbed free energy profile of the molecule 



(6) 

I 

odic loading, the best results were obtained for the opti- 
mal driving frequency uj* ~ 1.2 x 10 3 1/s, which yields 
an error of a 2 ~ 0.9 (k B T) 2 . This frequency is some- 
what smaller than the spontaneous transition rate under 
preloading, which is about 7 x 10 3 1/s for our model 
system. By analyzing the work distribution we have 
convinced ourselves that uj* indeed corresponds to non- 
equilibrium conditions. For both smaller and larger fre- 
quencies than uj* the quality of the reconstructed data 
becomes worse as expected from our reasoning above. 

For linear loading, the error increases for increasing 
driving velocity v, as expected. Since we fixed the total 
measuring time and the total number of trajectories, the 
least possible velocity for overcoming the barrier is 
v ~ 100 nm/s. The smallest error of a 2 ~ 1.7 (k B T) 2 , 
however, was observed at a larger, optimal velocity v* ~ 
200 nm/s, with an error bar of 0.03. For the unbiased 



G(z) = -B- 1 ] y {6(z~at))eM-PW(t))) /sr cxp^f3V (z,x(t))} 
{ ' ' V (exp(-/W(i))) / 4- (exp(-f3W(t))} 
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stants on each part, similarly as shown in [28( 
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Figure 3: Mean square error a and error bars for recon- 
structed free energy profiles by using the linear ramp J2) as 
a function of v (top scale, +), and the periodic ramp © as 
a function of u) (bottom scale, 0) (compare Fig.|2J. Both 
methods use 10 runs with 10 trajectories each. The optimal 
velocity v* and frequency u*, respectively, indicated by the 
vertical dashed line, correspond to the smallest error. 

comparison in Fig.|2 we have chosen the optimal values 
lo* and v* from the data shown in Fig.|3 

Discussion and summary - We have proposed a new 
method for recovering the free energy profile G(z) of 
biomolecules in single molecule pulling experiments by 
combining the new periodic ramp Q with Jarzynski's 
identity for recovering equilibrium properties from non- 
equilibrium data. The simulated data in Fig.[3show that 
the periodic ramp delivers significantly more accurate 
data than the traditional linear loading (J2J, under the 
constraint of equal experimental effort. An additional 
advantage of the periodic ramp is the fact that the mea- 
suring time may be chosen as long as necessary to col- 
lect the sufficient amount of data to recover the barrier 
regions of the free energy profile, which are hard to de- 
termine by previous methods. 

The driving frequency u> and preloading offset x$ of 
the periodic ramp J3J provide handles to optimize its 
performance (see Fig.OJ. As our model case study has 
shown, the frequency should be large enough to drive 
the system out of equilibrium, but still smaller than the 
spontaneous transition rate under preloading. Whether 
the optimal frequency is always of this order remains to 
be seen in more comprehensive systematic investigations 
left for future work. 

Our method can be extended to probing free energy 
profiles with different transitions between a number of 
metastable states. If these states are sufficiently sepa- 
rated, by several nanometers and energy barriers of sev- 
eral fcgT, each transition can be selected by using suit- 
able values for preloading offset xq and amplitude a in 
© ■ Subsequently, these parts of the free energy profile 
may be fitted together by adjusting free additive con- 
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